subroutine vmd_init_ex(kXvec, kZvec, pvec, vvec, u, v0, &             ! inputs
                       ExFilter, ExpKVdt, ExpKz, ExpKv, ExpKx)        ! outputs

  use dimensions, only : jx, jz, jv, jDirection,          &   ! domain size in r-space
                         nkx, nkz,                        &   ! limits      in k-space
                         dv                                   ! resolution  in r-space


  use parameters, only : phi0, kz, vPhase,  Dxxxx, Dvvvv, &            ! real,     phys
                         ny,  iWaveFrame,                 &            ! integer,  phys
                         Vmax, dt, KxFilter, &                         ! real,     num 
                         pi, iunit, isign                              ! constants

  implicit none

  real, intent(in)      :: kXvec(jx), kZvec(jz), pvec(jv)
  real, intent(in)      :: vvec(jv)
  real, intent(in)      :: u, v0

  real,    intent(out)  :: ExFilter(jx)
  complex, intent(out)  :: ExpKVdt
  complex, intent(out)  :: ExpKz(jz)
  complex, intent(out)  :: ExpKv(jz,jv)
  complex, intent(out)  :: ExpKx(jx,jDirection)


!-- local variables ---

  integer              :: ix, iz, iv, iDirection  ! counters in r-space
  integer              :: ikx, ikz, ip            ! counters in k-space
  real                 :: uVec(-1:1)              ! (-u, 0, u)
  real                 :: capOmega


!-- initialization of --------------------------------------------------------------------------------

  uVec(-1) = -u; uVec(1) = u ; uVec(0) = 0

  ExFilter=1.0
  do ikx=1,jx
     if( abs(kXvec(ikx)) > KxFilter ) ExFilter(ikx) = 0
  enddo

  ExpKz(:) = cexp( iunit * kZvec(:) * V0 * dt )  !  complex ExpKz(jz)
      

  do ikz=1,nkz
     ExpKv(ikz,:) = cexp( -iunit*vvec(:)*kZvec(ikz)*dt - Dxxxx * kZvec(ikz)**4 *dt    )
  enddo
      

  ! in lab frame, E0z has angular frequency vPhase*kz
  ExpKVdt = 1
  capOmega = vPhase*kz
  if(  iWaveFrame==0 ) ExpKVdt = exp( - iunit*0.5*dt*capOmega  )
                                                                            

  !  data iSign(1), iSign(2)/ 1, -1 /
  do iDirection = 1,2
     ExpKx(:,iDirection) = exp( - iunit * kXvec(:)*u*iSign(iDirection) * dt  )
  enddo


!-- Legacy output ----------------------------------------------------------------------

      open(5, file='output', position='append')
     
      write(5,*) "phi0=amplitude of external potential=", phi0, &
               "  kz/kDebye=", kz, &
               "  phase velocity/ve=", vPhase, &
               " dkz=kz /", ny, &
               "  jz= ", jz, "     ny= ", ny, "    jv= ",jv

      write(5,*) "Vmax/ve= ",Vmax,"  Omegape*dt= ",dt
          
      write(5,*)"number of Fourier modes of a given sign=",-1+nkz/3
      write(5,*)"dt figure of merit for advection = dt kmax Vmax =", &
                 dt*( (nkz/3)-1 )*kz*max(abs(vvec(1)),abs(vvec(jv)))     
      write(5,*) ' Dxxxx * k_zmax^4 / (k_zmax * dv) =', &
                 Dxxxx * (kZvec(nkz/2  - nkz/6 -1))**3 / dv, &
                 '  should be > 1 to resolve 1/kv  resonance'
       
      write(5,*)
      write(5,*)'Kx filter if |kx| > ', KxFilter ; 
      write(5,*) 'max |k| that survives de-Aliasing =',  kxVec(nkx/2  - nkx/6 -1)
      write(5,*)
      write(5,*)'kx '
      write(5,*) kxvec(:)
      write(5,*)
      write(5,*)'ExFilter '
      write(5,*) ExFilter(:)
      write(5,*)
      
      write(5,*)
      write(5,*)'kXvec(:) * dt'
      write(5,*) kXvec*dt
      write(5,*)
      write(5,*)
      write(5,*)'[kXvec/kz]^2' ; write(5,*) ( (kXvec(ix)/kz)**2, ix=1,jx/2 )
      write(5,*)
      write(5,*)
      write(5,*)'[kXvec/kz]' ; write(5,*) ( (kXvec(ix)/kz), ix=1,jx/2 )
      write(5,*)
      
      
      write(5,*);
      write(5,*) "Dvvvv x pvec(max)^4 x dt=",  Dvvvv * pvec(jv/2)**4 * dt, &
                 "   Dvvvv x pvec(max)^4=",    Dvvvv * pvec(jv/2)**4
      write(5,*) "Dxxxx x kZvec(max)^4 x dt=", Dxxxx * kZvec(nkz/2)**4 * dt,&
                 "    Dxxxx x kZvec(max)^4=",  Dxxxx * kZvec(nkz/2)**4
      write(5,*) "kz(max)*v(max)=",  kZvec(nkz/2)*max(abs(vvec(1)),abs(vvec(jv))), &
                 "    p(max)=", pvec(jv/2)
      write(5,*)

      close(5)
!-------------------------------------------------------------------------------------

end subroutine vmd_init_ex
